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STEADILY TRANSLATING PARABOLIC DISSOLUTION FINGERS 

PAWEL KONDRATIUKt AND PIOTR SZYMCZARt 

Abstract. Dissolution fingers (or wormholes) are formed during the dissolution of a porous rock 
as a result of nonlinear feedbacks between the flow, transport and chemical reactions at pore surfaces. 
We analyze the shapes and growth velocities of such fingers within the thin-front approximation, in 
which the reaction is assumed to take place instantaneously with the reactants fully consumed at the 
dissolution front. We concentrate on the case when the main flow is driven by the constant pressure 
gradient far from the finger, and the permeability contrast between the inside and the outside of the 
finger is finite. Using Ivantsov ansatz and conformal transformations we find the family of steadily 
translating fingers characterized by a parabolic shape. We derive the reactant concentration field 
and the pressure field inside and outside of the fingers and show that the flow within them is uniform. 
The advancement velocity of the Anger is shown to be inversely proportional to its radius of curvature 
in the small Peclet number limit and constant for large Peclet numbers. 

Key words, free boundary problems, dissolution, convection-diffusion-reaction, porous media 

AMS subject classifications. 35R35, 86A60, 76S05, 35Q35, 35Q86, 30E25 

1. Introduction. Chemical erosion of a porous medium by a reactive fluid is 
relevant for many natural and industrial applications and is a crucial component 
in a number of geological pattern formation processes [36, 29]. Networks of caves 
and sinkholes are formed by the dissolution of limestone by C02-enriched water in 
karst areas [40, 50], the ascending magma dissolves the peridotite rocks, leading to 
formation of porous channels [2], matrix acidizing is used by petroleum engineers to 
enlarge the natural pores of the reservoirs [45] - in all of these processes the interplay 
of flow and reaction in an evolving geometry results in spontaneous formation of 
intricate patterns. Reactive flow systems with strong, nonlinear coupling between the 
transport and geometry evolution have been investigated in a number of experimental 
[15, 16, 26] and theoretical [42, 20, 18, 32] studies, the latter often combined with 
numerical simulations. 

A locally increased dissolution rate can make the rock more porous and thus 
more permeable. This increases the flow and reactant transport, further enhancing 
the local porosity increase. Such a positive feedback eventually leads to the creation 
of highly localized flow paths, which go by different names depending on the field. In 
petroleum industry, they are dubbed “wormholes” [23] due to the resemblance with 
the tunnels dug by worms, geologists call them “solution pipes”, “karst funnels” or 
“geological organs” [30, 53, 17], whereas pedologists - “soil tongues” [54]. Examples of 
such structures exposed in the limestone quarry in Smerdyna, Poland, are presented 
in Fig. 1.1. 

The initial stages of the wormhole development due to the breakup of the ini¬ 
tially planar reaction front (referred to in the literature as the “reactive-infiltration 
instability”) have been extensively studied [8, 39, 48, 9, 22, 52] and are now quite 
well understood. Much less is known, however, about the nonlinear regime, when the 
initial perturbations of the interface are transformed into finger-like structures that 
advance into the system [15, 21, 41, 11, 49]. An important question in this context 
concerns the existence of a self-preserving form that might model the growing tip. 
Such invariantly propagating forms have been found in other pattern forming sys- 
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Fig. 1.1. Solution pipes in the limestone quarry in Smerdyna (Poland), formed in limestone 
bedrock. They are postulated to have been formed by a dissolving action of meltwater during Elsterian 
deglaciation [37, 53]. The brown rims in the clay filling the pipes are formed as a result of clay and 
iron oxide accumulation due to the illuviation processes [47j. The clay particles are transferred by 
water from the upper parts of the soil and then flocculated at the clay-limestone boundary, where 
the pH changes from mildly acidic to alkaline. 


terns, e.g. the SafFman-Taylor finger in viscous fingering or Ivantsov paraboloid in the 
dendritic growth [46, 24, 28, 1, 14]. Our goal in this paper is to find such solutions 
for a dissolving porous medium. 

In the context of dissolution, the propagation of individual fingers was studied 
in the petroleum engineering [27, 6, 11]. In these works, simple geometric models 
of wormhole shapes were adopted, which are not preserved during the growth. The 
closest in spirit to the present study is the work by Nilson and Griffiths [38]. They 
tackle the problem of steadily propagating dissolution forms and find them to be 
parabolic (in 2d) or paraboloidal (in 3d). However, there are significant differences 
between their approach and the present one. We discuss those in more detail at in 
Sec. 5. Here we will just note that in [38] the reactant concentration field is not 
resolved, instead the dissolution front velocity is taken to be proportional to the local 
fluid velocity. 

Additionally, the pressure drop inside the finger is neglected, which is only justified 
if the permeability of the dissolved phase is much higher than the permeability of the 
primary rock. In that way, the problem is effectively been reduced to a one-phase 
problem, requiring the solution of the Laplace equation for pressure in the region 
outside the finger only. As discussed in [25, 12] such “one-phase” growth problems 
are in general much easier to solve than the “two-phase” problems in which the flow 
fields both inside and outside the finger need to be found, which is also the case in 
the present study. On top of it, to find the finger advancement velocity, we need to 
solve not only for the pressure, but also for the reactant concentration field, which 
makes the task at hand more challenging. 

The paper is organized as follows. In §2 the general equations governing the 
dynamics of the dissolving porous rock are briefly recalled and then put in the di¬ 
mensionless form in §3. The core of the paper are sections 4 and 6, where - after 
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Fig. 2.1. Geometry of the system: a reactive fluid is injected from the left side and dissolves the 
porous matrix through chemical reactions. Once the dissolution is complete, the reaction front, shown 
by the solid line, advances into the matrix (dashed line), separating the dissolved, upstream domain 
(O) of porosity (fi and the undissolved, downstream domain (Cl) of porosity cj>g, complementary to 

n. 

the application of the Ivantsov ansatz - we obtain the family of steadily translating 
dissolution fingers and show that they are parabolic/paraboloidal. Let us reiterate 
that despite a formal similarity to the Ivantsov forms, the physics of growth is quite 
different here, as it involves an interplay between the flow field and solute concen¬ 
tration field in contrast to the dendritic growth which is controlled by a single field 
(temperature only). 

2. The model of matrix dissolution. When a porous matrix is infiltrated by 
an incoming flux of reactive fluid, a front develops once all the soluble material at 
the inlet has been dissolved. This front propagates into the matrix as illustrated in 
Fig. 2.1. Upstream of the front, all the soluble material has dissolved and the porosity 
is constant, (p = (f>i. Ahead of the front the porosity decays gradually to its value in 
the undissolved matrix, (p = (po- The front is initially planar but eventually breaks 
up because of a positive feedback between flow and dissolution, which amplihes any 
small variation in the porosity field [8]. 

Let us briefly recall the equations for the dissolution of a porous matrix. Rate 
of groundwater flow through the porous medium is taken to be proportional to the 
pressure gradient (Darcy’s law). 


V = 



( 2 . 1 ) 


where v is the Darcy velocity, p is the porosity and K(p) is the permeability. We also 
assume that the Darcy velocity field is incompressible, 

V • V = 0, (2.2) 


neglecting contributions to the fluid volume from reactants or dissolved products. 
Under typical geophysical conditions, dissolution is slow in comparison to flow and 
transport processes; we can therefore assume a steady state in both the flow and 
transport equations. 

The transport of reactants is described by a convection-diffusion-reaction equa¬ 
tion, 


V • (vc) -y -{DP- Vc) = -R, 


(2.3) 
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where c denotes the concentration of the reactant and R{c) is the reactive flux into 
the matrix. In the upstream region, where all the soluble material has dissolved {(j) = 
(pi), the reaction term vanishes and the transport equation reduces to a convection- 
diffusion equation. 

In the derivations below, we adopt a thin-reaction-front approximation, in which 
the transition zone over which the porosity changes is assumed to be infinitely thin 
[ 8 , 39]. As shown in [51], this approximation is applicable whenever the reactant 
penetration length is short in comparison with the diffusive length scale D/v, char¬ 
acterizing the decay of the concentration in the upstream region. In practice, this 
happens whenever the parameter H = Dk/v^ ^ 1 [51], where k is the dissolution 
rate (assuming the 1 st order reaction kinetics); thus the reaction rate needs to be 
sufficiently high and/or the flow rate — sufficiently low. As a consequence, we obtain 
a Stefan-like problem in which the space is divided into two domains: the dissolved, 
upstream domain (If) of porosity pi th® undissolved, downstream domain (fl) 
of porosity po, complementary to fl. The reaction front (cAI) advances with velocity 
proportional to the flux of the reactant at a given point 

Un = -—D{Vc)^ (2.4) 

Cin 

where subscript n represents the component normal to the interface, is the inlet 
concentration of the reactant and the acid capacity number 7 = Cin/vCsoi[Pi — Po) 
is defined as the volume of rock (of molar concentration Csoi) that is completely 
dissolved by a unit volume of reactant (of molar concentration Ci„). Finally, v is the 
stoichiometric coefficient in the dissolution reaction (number of moles of the reactant 
necessary to dissolve one mole of the rock). 

The flux of reactant in (2.4) is taken to be purely diffusive, since the thin-front 
approximation implies vanishing of c at the interface, as it is fully consumed there: 

c|9n(t) = 0- (2.5) 

Because in both domains the porosity (and hence permeability) is constant, 
Eqs. (2.1-2.2) can be combined to yield the Laplace equation for the pressure 

= 0 (2.6) 

in both n and II. 

The above equations are supplemented by boundary conditions on the velocity 
and concentration field as x —>■ 00 : 

v(x —>■ 00 ) = voex, c(x 00 ) = 0, (2-7) 

which accounts for the fact that far from the dissolution front the flow becomes 
uniform and the reactant concentration vanishes. The conditions at x —>■ —00 are 
somewhat more subtle. If the reaction front assumes such a form that sufficiently far 
upstream everything is dissolved, i.e. 3xo (Vx < Xq, {x,y) € H), then it is possible to 
impose 


dxVx{x - 00 ) = 0 


( 2 . 8 ) 


on the flow field and 


c(x -)► - 00 ) = dn, 


(2.9) 
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on the concentration field; the former representing the condition that the flow becomes 
uniform far upstream, the latter corresponding to the reactant concentration imposed 
at the inlet. 

However, if the undissolved phase (fl) extends towards x —>■ —oo (as it is the 
case for a solitary dissolution finger surrounded by an undissolved matrix), then the 
condition analogous to (2.9) should be imposed only within the finger, and not on the 
entire a; —>■ — oo line. We will come back to this issue later on. 

Additionally, both the pressure and the normal component of the Darcy velocity 
need to be continuous across the reaction front dfl: 

p\dn{t)-= p\dn{t)+, (2-10) 

i^n|an(t)-= Tn|on(t)+ (2-11) 

3. Scaling of the variables and the limiting cases. The dissolution equa¬ 
tions (2.1-2.3) can be simplified by scaling the velocity, concentration, and porosity 
fields by their characteristic values 


v = v/uo, c = c/c™, (f>=4- —(3.1) 

01 — 00 

where the scaled, dimensionless variables are marked by hats. Additionally, we scale 
the spatial coordinates by some length I, characterizing the wormhole, and time by 
r = l/jvo- Using I, we can scale the pressure as follows 


^ KiM 
p = - rP- 

Vopl 

With these scalings, the governing equations take the form 


in the downstream one. In the above, 

Pe = 


VqI 


0iD(0i)’ 


(3.2) 


= 0 

r S n{i) 

(3.3) 

AtVp • Vc -1- Pe = 0 

r £ Q{t) 

(3.4) 

in the upstream region and 

<l> 

K5 

II 

o 

r £ n{i) 

(3.5) 

o 

II 

r £ Cl{i) 

(3.6) 


(3.7) 


is the Peclet number which measures the relative magnitude of convective and diffusive 
effects on the length scale I, whereas 


_ A^(0i) 

K{M' 


(3.8) 


is the ratio of permeability between the domains. In typical geological conditions 
the permeability is an increasing function of porosity, thus k > 1. The boundary 
conditions (2.7-2.9) in the scaled variables take the form 


v(x oo) = e, 


c{x —>■ oo) = 0, 


(3.9) 
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and 


{dxVx)(x —>■ —oo) = 0, c{x —>• —oo) = 1, (3.10) 

which, as before, need to be supplemented by the continuity conditions for the pressure 
p and the normal component of the velocity Vn across the interface^ at the boundary 
(reaction front) Finally, the condition (2.4) for the front advancement velocity 
takes the form 

17„= Pe-^(Vc)„ (3.11) 


Noting that at the reaction front the concentration satisfies 

A ■' 3c 

= (3.12) 

ot 

one can rewrite (3.11) in terms of c only: 

5c(r) ^ Pe~^|Vc(r)p = 0 r S 50(t) (3.13) 

ot 

It is relatively straightforward to derive the one-dimensional solutions of Eqs. (3.3- 
3.6), corresponding to the planar reactive front propagating with a constant velocity, 
Uq. Assuming that both c and p are the functions of x only, leads to 

..... f —K~^X X <0 /n 1 

pM = I i <=>■“> 


where x = x — Uot is the coordinate moving with the front. For the concentration we 
get an exponentially decaying, diffusive profile 


c(i) 


l-eP®“ i<0 
0 i > 0. 


(3.15) 


Finally, inserting the above into the front velocity condition (3.11) we obtain the 
result Uq = 1, i.e. in the units chosen the planar reaction front advances with a unit 
velocity along the flow direction. 

4. Stationary dissolution fingers in 2D. In this section, we derive the shape 
of a single dissolution finger propagating invariantly towards the undissolved matrix. 
We begin with the two-dimensional case. 

At the boundary between the phases, Eq. (3.13) is satisfied. Without any loss in 
generality, one might extend this equation to the whole domain 17 as 

Br 

F(f,t)^-|Vcp = 0, (4.1) 

at 

with F formally given as F = |Vcp/i9jC and satisfying = Pe. To proceed, 

we use the ansatz due to Ivantsov [28, 1] by assuming that the unknown F function 
dependence on the spatial and time coordinates is of the form 

F(f,t) = F(c(f,7)), (4.2) 


^Note that v = —AcVp on Q and v = —Vp on Q. 
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so that 


F(c)g-|Vcp = 0 (4.3) 

at 

is satisfied in the whole domain D,. Obviously, F must satisfy 

F(c = 0) = Pe (4.4) 

for the consistency with the boundary condition (3.13). Using this idea, Ivantsov has 
found his famous set of steady-state solutions to the dendrite growth problem in a 
supercooled melt. Note that assuming the ansatz (4.2) might lead to losing some of 
more general solutions of the problem. 

Let us now suppose that our problem has a stationary solution: a dissolution 
finger moving invariantly in the x direction with velocity U. Rewriting Eq. (4.3) in 
the moving coordinate system gives 


-UF{c)—-\Vc\^ = 0 


(4.5) 


or 


Vc • (Vc-b UF(c)Vi) = 0. (4.6) 

Let us note that both the Laplace equation (3.3) and the advection-diffusion equa¬ 
tion with potential flow (3.4) are conformally invariant. A large class of conformally 
invariant, non-Laplacian physical problems has been identified by Bazant in [3, 4] 
and involves such phenomena as nonlinear diffusion or advection or electromigration 
coupled to diffusion. The conformal invariance of these processes provides an effec¬ 
tive way of solving these problems [10]. In the context of growth processes, these 
techniques have been used to track the evolution of the interface in solidification and 
melting under the action of a potential flow [34, 19, 33, 14, 13, 3, 5]. However in these 
works the growth has been taking place in external flows (with the fluid outside of the 
growing object), in contrast to the present case where the medium is porous, and the 
driving flow goes through the growing finger and then into the undissolved medium. 

In our case, Eq. (4.6) is also written in a conformally invariant form. This allows 
us to make a conformal coordinate transformation {x,y) —> (^,??) such that c = c(^) 
(i.e., the isolines of c coincide with the curves ^{x,y) = const). Since Eq. (4.6) is 
conformally invariant, we get 


The immediate conclusion is that ^ must be a function of ^ only^, hence 


(4.7) 


x{^,r]) = $i(^) -b 4>2(??)- 


(4.9) 


^Note that if one sought a stationary solution without having assumed the Ivantsov ansatz, one 
would obtain at this point 


dc 


^'T 

+ UF{^,v)—=0 

9 ^ 


(4.8) 


which does not provide any information about the mapping. 
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The Cauchy-Riemann conditions, 


dy dx dy dx 
drj drj' 


(4.10) 


lead to the Laplace equation for x, = 0, the general solution of which, compatible 
with (4.9), is 


v) = ((C - -{v- V*f) ■ 


(4.11a) 


Using again (4.10) allows us to find y{^,ri) in the form 

ii) = p(C V*) + y*- 


(4.11b) 


Without the loss of generality we can set = ? 7 , = j/* = 0. Furthermore, we can 
assume that the parabola characterized by ^ = 1 corresponds to the dissolution front 
(any other parabola can be scaled onto it by an appropriate choice of parameter p). 
Note that the parameter p corresponds to the (dimensionless) radius of curvature 
of the parabolic finger. This suggests that the (dimensional) radius of curvature, p, 
could be used as a length parameter in our problem. Taking I = p leads to p = 1 
which simplifies the subsequent analysis. 

Our final conclusion is the following: the only conformal mapping {x,y) —^ {CtV) 
which makes the concentration field in the dissolution finger a function of one variable 
only, c{x,y) = c{^{x,y)), is (up to translation)^ 




(4.12) 

(4.13) 




which is the usual transformation to the parabolic coordinates. 

It is worth noting that parabolic geometries have also been obtained by other au¬ 
thors studying steadily translating growth forms. The classical example are Ivantsov’s 
dendrites [28]. In fact, Adda Bedia and Ben Amar [1] have shown that, at least in 
the context of steady state dendritic growth at zero surface tension, the parabolic 
solutions are the only admissible solutions once the Ivantsov ansatz is adopted and 
that even when one modifies this approach by making a more general ansatz, no new 
forms are found. Although the present case is of a more complicated nature, it is 
possible that also here the adoption of Ivantsov ansatz necessitates the appearance of 
the parabolic forms. 

Once the shape of the finger is known, we can find the concentration and pressure 
both inside and outside of it. 

4.1. Solution inside the finger. Inside the dissolution finger (i.e., for ■C < 1) 



(4.14) 


d'^p d'^p 
drf 


+ TTW — 0 


(4.15) 



(4.16) 


^In fact, we also assume the Ivantsov ansatz (Eq. (4.3)). 
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Fig. 4.1. Schematic representation of the boundary conditions for the concentration field c 
inside the parabolic dissolution finger. The finger is represented in the Cartesian coordinate system, 
with the X axis pointing rightwards and the lengths scaled by the radius of curvature p. The outer 
parabola, denoted by ^ = 1, constitutes the boundary of the finger and is characterized by vanishing 
of the reactants concentration, c{^ = 1) = 0. At the inner parabola of ^ = Co (here Co = 0-lJ the 
boundary condition (2.9) is imposed, c(C = Co) = 1- Additionally, the line C = 0 is shown. 


Additionally, at the front itself, the reactant concentration vanishes, i.e. c(^ = 
1) = 0. It is harder to formulate the second boundary condition for the concentration, 
necessary to close the equation. As mentioned in §2, the upstream condition in the 
form of Eq. (2.9) is inconsistent with the presence of the finger extending towards 
X -A — oo. Instead, we need to prescribe the concentration within the finger. In 
accordance with the assumption that c is a function of ^ only, we impose c = 1 on 
one of the parabolas, dehned by ^ (see Fig. 4.1). 

To proceed, we conclude from Eq. (4.14) that d^p must be a function of ^ only. 
Along with the fact that p satisfies the Laplace equation (4.15) this gives, similarly 
as before, 


=-^A(C^ - ?7^) + Ai^ +A 277 (4.17) 

where the constant term has been dropped. However, the constants Ai and A 2 have 
to be zero, otherwise the flow v = —kVp would diverge at f = 0. Inserting Eq. (4.17) 
back into (4.14) yields for the concentration field 


. erfi(a^) — erH(Q;) 
erfi(a^o) ~ erfi(a) ’ 


(4.18) 


where a = Pe and erh(a;) = e‘ dt. 

Finally, Eq. (4.16) allows us to hnd the velocity of the propagation of the finger 


U = - 


c'iO 


F{c{^))d^x 


2a e° 




n( Pe (erfi(a) — erh(afo)) ’ 


(4.19) 


where we have used the fact that E(c(^ = 1)) = Pe. 

4.2. Solution outside the finger. Outside the dissolution finger, in the undis¬ 
solved matrix, the pressure still obeys the Laplace equation 


d'^p d'^p 

9^2 Qfj2 


(4.20) 
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The continuity conditions for p and the yield here 




- di 


?^i+ 


This is supplemented by the boundary condition (2.7) 


—Vp 




- - V X , 


which can be written in the (^, 77 ) coordinates as: 


y ^ dp 
5—>oo ^ 


lim = 1. 

5^00 77 drj 


The solution of Eq. (4.20) fulfilling (4.21-4.24) is 

- («:- 1)(^- !)■ 


(4.21) 

(4.22) 

(4.23) 

(4.24) 

(4.25) 


Additionally, the condition (4.22) gives A = 1 for the constant appearing in the 
solution of the internal problem (Eqs. 4.17-4.18). 

4.3. Summary. The pressure field in the system is 


P(C,»7) 


-2(^^-77^) forO<^<l 

-l(C^-?7^)-(«-l)(^-l) for^>l. 


(4.26) 


Fig. 4.2 shows the pressure fields, plotted for various permeability contrasts k, in the 
physical (a;, y) plane. Note that inside the finger the isobars are parallel to each other 
and oriented along the y direction, which is a manifestation of the fact that the flow 
is uniform there: 


P{x,y) = -x, refl, (4.27) 

'v{x,y) = KBx, reft. (4.28) 

In the dimensional variables, the (constant) Darcy velocity inside the finger is 

Vin = KVoBx (4.29) 

Outside the finger, in the undissolved domain D, the pressure field is 

p{x, y) = —X — (k — 1) + X — 1^ . (4.30) 

The term linear in x dominates far from the boundary i9D, which agrees with the 

far-downstream condition for the flow field (2.7). Note that both the pressure field 
inside the finger (4.28) as well as that outside (4.30) obey the upstream boundary 
condition ( 2 . 8 ). 

The above-obtained pressure field is analogous to that derived by Kacimov and 
Obnosov [31], who have been analyzing the groundwater flow in a porous medium 
with a parabolic inclusion. In fact the pressure problem solved in [31] is more general 
than that considered here, since it involved the far-held how velocity oriented at an 
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Fig. 4.2. Pressure field inside and outside of the two-dimensional stationary dissolution finger 
(marked by the thick parabola), for different values of permeability contrasts ft, ft = 2 (a), k = 10 (b), 
and ft = 100 (c). The (dimensionless) pressure drop between the neighboring isobars is 2 (a), 5 (b), 
and 50 (c). 


arbitrary angle with respect to the parabola axis, not necessarily parallel as in (2.7). 
Interestingly, in each case the flow within the parabola is found to be uniform. This is 
the parallel of the classical result by Poisson and Maxwell who noted that the electric 
field inside elliptic or ellipsoidal inclusions is uniform [44, 35], see also Carslaw and 
Jaeger [7] for the interpretation of this fact in the context of heat transfer. 

Finally, let us note that the higher the permeability contrast k, the more disturbed 
is the pressure field by the presence of the finger and the more aligned the isobars 
along the reaction front. In the limit of k —^ oo, the dissolution front dil would 
corresponds to one of the isobars. 

The concentration c is nonzero only inside the dissolution finger and depends on 
the ^ coordinate in the following way. 


5(0 = 


erfi 


(\/- erfi 




erfi 


(v^I^PeO) 


— erfi 




for 0 < C < 1 


(4.31) 


where, as mentioned above, the inlet concentration, c = 1, has been imposed on the 
parabola ^ = io (0 < < !)■ As a special case, one can choose .Jo = 0, which yields 


c(J) = l- 


erfi iftPe 
erfi (^iftPe) 


for 0 < J < 1. 


(4.32) 


The inspection of Eq. (4.31) reveals that the relevant parameter controlling the con¬ 
centration field is the product 

Pe' = KPe=-^. (4.33) 

01U 

with Vin given by Eq. (4.29). This parameter is interpreted as a Peclet number related 
to the flow within the finger. 

The impact of this parameter on the concentration profiles is analyzed in Fig. 
4.3. For small values of Pe^ Eq. (4.31) is approximated by 


c 


e-1 

Jo-1 


for Jo < J < 1- 


(4.34) 
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Fig. 4.3. Reactant concentration field c inside the two-dimensional stationary dissolution finger 
(marked by the thick parabola), for Pe' = 0.1 (a), Pe' = 5 (b), and Pe' = 10 (c). The boundary 
condition for the concentration (2.9), has been imposed on the line ^ = 0. The (dimensionless) 

concentration drop between the neighboring isolines of c (thin parabolas) is set to 0.1. The boundary 
of the finger is also an isoline of the concentration (for c = 0). 


i.e. the concentration depends linearly on the parabolic coordinate The gradient of 
the concentration field is then relatively uniform (cf. Fig. 4.3a). This situation can be 
referred to as the “diffusive” regime. In the opposite case of large Pe' (the “convective” 
regime), the concentration boundary layer is formed, with uniform concentration c = 1 
in the body of the finger and high gradients near the boundary (Fig. 4.3c). In this 
limit the solutions found by us acquire direct physical relevance since a somewhat 
artificial choice of becomes irrelevant here, as long as the parabola C = Co lies in 
the bulk. The characteristic width of the boundary layer is (^Pe') , therefore the 

condition 


Co "C 1 



(4.35) 


guarantees that Co lies in the bulk and thus its precise value does not affect the 
concentration profiles. This is confirmed by the comparison of c profiles for Pe' = 10 
and Co = 0.05 with those for C = 0.2, which are found to be almost indistinguishable, 
as shown in Fig. 4.4. 

Finally, the advancement velocity of the dissolution finger is 


£'(0 


•\/2k exp 

(Fe') 


F{c{C))d^x 

Vtt Pe 

erfi 


- erfi 

(v^fo)) 


(4.36) 


with the two limiting cases (corresponding to small/high values of Pe') given by 


lim Pef7 = (1-Co)"\ (4.37) 

Pe'->0 

lim U = K. (4.38) 

Pe'— 


The velocity scale is I /t = ■jvo. Thus in the limit of small Pe', the dimensional finger 
advancement speed approaches 


'yvo _ l4>iD 

Pe(l-Co) " P(l-Co)’ 


U 


(4.39) 
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Fig. 4.4. (a,b) Reactant concentration field c inside the two-dimensional stationary dissolution 
finger in the “convective” regime, with Pe' = 10. The thick solid lines mark the finger boundary and 
the C = ^0 parabola with = 0.05 (a) and = 0.2 (b), on which the inlet concentration has been 
imposed, c(^ = Co) = 1- The (dimensionless) concentration drop between the neighboring isolines of 
c (thin parabolas) is set to 0.1. The dashed line marks the parabola tj = 2 (locally orthogonal to the 
isolines of concentration) parametrized by the arc length s. The reactant concentration profile along 
this parabola, c(s), is plotted in panels (c) and (d). 



i.e. it is inversely proportional to the radius of curvature of the tip. Such a behavior 
is in full analogy to the Ivantsov result for the dendrite propagation speed [28]. On 
the contrary, for large Peclet numbers the propagation velocity becomes constant with 
respect to the curvature of the finger tip and simply proportional to the Darcy velocity 
inside the finger 


U Ri fiVoK = IVin- (4.40) 

Notice the inherent indeterminacy of the general result (4.36): we have obtained the 
relation between the tip radius and the propagation speed with no means of evaluating 
each of these quantities independently. This is connected with the scale-invariance 
of the original equations (2.1-2.4) and can be lifted only after the introduction of an 
additional lengthscale. Again, this feature is shared by both the Ivantsov parabolas 
and the Saffman-Taylor fingers, which also represent families of growth forms and 
need the short-scale regularization mechanism such as surface tension or kinetic un¬ 
dercooling for the selection of a particular tip radius and advancement velocity. In 
the case of dissolution patterns, the additional lengthscale might originate from the 
reaction front width (assumed to be zero in our analysis) or from the interactions 
with other dissolution fingers: in real system, the finger is never infinite, but always 
constrained by the presence of its neighbors. 

5. Comparison with the Nilson & Griffiths model. At this point it is 
worth noting the main differences between our finger growth model and the one of 
Nilson and Griffiths [38]. First of all, we differ in the boundary conditions at infinity. 
Whereas we assume that the entire system is flushed with a reactive fluid (with 
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constant pressure gradient both at a; —>• —oo and x —>■ oo), Nilson & Griffiths assume 
constant pressure at a; —> oo, which means that the fluid far from the tip of the finger is 
kept immobile. Since at the same time the reactive fluid is constantly injected into the 
system, displacing the original pore fluid, the latter becomes compressed, which does 
not seem to be a realistic assumption under typical groundwater pressures. In another 
case considered in Ref. [38] the medium is supposed to be initially unsaturated, with 
no fluid present. Then the reactive fluid is entering the system, creating a finger 
and flooding the medium in front of the finger tip, forming a ’product layer’. Due 
to the absence of imposed pressure gradient, the fluid in the product layer gradually 
slows down, as the thickness of the layer is increased until it reaches a stationary 
situation in which both finger and the product layer advance with equal speeds into 
the medium. Again, this is a fundamentally different situation from the one considered 
in the present work, where the entire medium is flushed with a reactive fluid due to 
the imposed pressure gradient. In particular, the solution tubes of Fig. 1.1 are formed 
under such conditions, with the gravity imposing the pressure gradient across the 
entire system. 

There is yet another difference between our approaches, this time related to the 
constraints on the values of the physical parameters, under which the finger model is 
supposed to work. Similarly to us, Nilson & Griffiths adopt the thin-front approxima¬ 
tion, stating that the reactant is entirely consumed as soon as it hits the finger bound¬ 
ary. As mentioned above, this corresponds to the assumption that H = Dk/vf^ » 1. 
At the same time, however, it is assumed in [38] that the concentration of the reac¬ 
tant inside the finger is constant, so that the dissolution front velocity is proportional 
to the local fluid velocity. This is only possible if the advection within the finger 
dominates over the diffusive transport, i.e. Pe^ = Vinl/D ^ 1. The double limit 
[Dk/vf^ » 1, Vinl/D ^ 1) puts rather stringent constraints on the admissible 
values of D, Vm and k — i.e. the flow rates should on one hand be large (to keep the 
transport advective within the finger), but on the other hand small enough to keep 
the reaction front thin. Conversely, in the present work, while keeping the front thin, 
we fully resolve the concentration field inside the finger, and thus obtain the finger 
propagation velocities over the entire range of Peclet numbers, c/. Eq. (4.36). At 
smaller Pe', the effects connected with the buildup of a diffusive layer of the depleted 
reactant appear [51], which have a significant impact on the wormhole propagation 
speed. 

The differences in the boundary conditions between our system and the one of 
Ref. [38] are reflected in different relations between the finger advancement velocity 
and its tip curvature. Nilson & Griffiths find that the finger velocity is inversely 
proportional to its radius of curvature. We get a similar result, but in the limit of 
small Pe' only, when the diffusive effects play a central role in the dynamics. These 
effects are absent altogether in the Nilson & Griffiths model. Conversely, we find that 
for large Pe^ the finger propagation velocity is independent of its curvature. 


6. Stationary dissolution fingers in 3D. In this section we study the geom¬ 
etry of three-dimensional dissolution fingers. Again, we assume the Ivantsov ansatz 
as well as the stationarity of the fingers, so the starting point of our investigations is 
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the following set of equations for pressure and concentration inside the finger 

kVp • Vc + Pe = 0 (6.1a) 

V2p=0 (6.1b) 

Vc • (F(c) t/ Vz + Vc) = 0 (6.1c) 


where we have assumed the flow along the z direction. As before, on the boundary 
of the finger (^ = 1) we put c = 0 and F = Pe. In the region outside the finger the 
pressure obeys the Laplace equation whereas the concentration vanishes. 

Let us switch to the cylindrical coordinate system, (x,y,z) —>■ (f,(p,z) and look 
for the solutions with rotational invariance, i.e. independent of cj). Eqs. (6.1) in the 
new variables take the form 


^ /V -'I dc 

KPeVrzP ■ '^rzC + -f T WT = 0 

r or 

r G n{i) 

(6.2a) 


r G fl{i) 

(6.2b) 

VrzC • (E(c) U VrzZ + ^rzC) = 0 

r G n(<) 

(6.2c) 

where Vrz = {dr,dz) acts as if the (r, z) were 

Cartesian coordinates. 

Similarly, 


. Analogously to the previous case we are looking for a conformal map 
(z,f) —such that c becomes a function of one variable only, c = c(C)- Again, 
from Eq. (6.2c) we draw the conclusion that such a map need to be a parabolic one, 

z = l{e-ri^) (6.3a) 

f = ^?7, (6.3b) 

with the radius of curvature of the finger chosen as the unit length. As before, we 
assume that the finger boundary coincides with the ^ = 1 paraboloid. 

Let us now use this mapping to solve the Laplace equation for the pressure in¬ 
side and outside of the dissolution finger. Since the Lame coefficients for the 

transformation (6.3) are equal to 



equation (6.2a) written in terms of ^ and rj becomes 

Pe'^c'{O+c"{O + jc'{O = 0- ( 6 . 5 ) 

Its form suggests that d^p is a function of ^ only and thus the pressure field p inside 
the dissolution finger is separable into 

+Y{r]). ( 6 . 6 ) 

The Laplace equation for the pressure inside the dissolution finger, Eq. (6.2b), in the 
rj coordinates transforms into 

d'^p d'^p \dp 1 dp 
9^2 drf ^ rj dr] 


0 , 0 <^< 1 . 


(6.7) 
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Fig. 6.1. Pressure field inside and outside of the three-dimensional stationary dissolution finger. 
The xz cross-section of the finger is shown, with the finger boundary marked by a thick parabola. 
The pressure field is plotted for different values of permeability contrast: k = 2 (a), k = 10 (b), and 
n = 100 (c). The (dimensionless) pressure drop between the neighboring isobars is 2 (a), 2.5 (b), 
and 10 (c). 


Inserting ( 6 . 6 ) into Eq. (6.7) gives the result that the pressure field inside the disso¬ 
lution finger should be of the form 

= -rf) + C^^^i + C2\mi + C^. 0 < ^ < 1 . ( 6 . 8 ) 

We set Cl = C 2 = 0, to avoid singularities at = 0 or r; = 0, while the constants Ai 
and Ca are to be determined from the boundary and the continuity conditions. Next, 
we look for the solution of the Laplace equation for the pressure field outside of the 
finger of the same functional form as that inside, 

P(C,? 7 ) =-^A2(C^-77^)-f Hiln^-l-DalnT?, ^>1. (6.9) 


where the constant term has been skipped. The continuity conditions and the far- 
downstream boundary condition for the flow are analogous to the two-dimensional case 
(Eqs. (2.7-2.11)). The downstream boundary condition yields A 2 = 1. The continuity 
of the pressure and the normal component of the flux at the finger boundary ^ = 1 
yields D 2 = C 3 = 0, Ai = 1 and Di = 1 — k. Eventually, the pressure field in the 
whole space is 




for 0 <^<l 

-l(?^-? 7 ^)-(K-l)lnC for^>l. 


( 6 . 10 ) 


This solution has been plotted in Eig. 6.1 for different values of the permeability 
contrast, k. The overall picture is similar to that in the two-dimensional case, with 
a uniform flow inside the finger and the flow disturbance around it relaxing to a 
uniform flow field far downstream. For small values of the permeability contrast k, 
the presence of the finger does not significantly disturb the isobars, whereas in the 
opposite case the isobars become almost parallel to the finger boundary. 

Using (6.10) one can solve Eq. (6.5) and find the concentration field inside the 
finger. Formally, the solution can be written in terms of the exponential integral 
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Fig. 6.2. Reactant concentration field c in the cross-section of the three-dimensional stationary 
dissolution finger, for Pe'=0.1 (a), 5 (h), and 20 (c). The parameter was chosen to be 0.1. The 
(dimensionless) concentration drop between the neighboring isolines is set to 0.1. 


function Ei^, 


e(e) 


Ei(lPe'a-Ei(lPeO 
Ei(iPe'Co^)-Ei(iPe')’ ^ 


( 6 . 11 ) 


where, as before, we set the condition c = 1 on the parabola C = Co (with 0 < Co < !)• 
The concentration fields for various parameters have been plotted in Fig. 6.2. Quali¬ 
tatively, the behavior is similar to the one observed for the 2D case. The “diffusive” 
regime is observed for small values of Pe^ (see Fig. 6.2a) and is characterized by 
relatively uniform distribution of the concentration gradient. In this case, the con¬ 
centration field can be approximated by® 


fi(C) 


InC 

In Co’ 


Co < C < 1- 


( 6 . 12 ) 


In the opposite case of high Pe' (see Fig. 6.2c), the boundary layer characterized by 
high gradients of the concentration c is formed near the finger boundary, while the in 
the interior of the finger the c field is almost uniform. The characteristic width of the 
boundary layer in the (C, ? 7 )-plane is, similarly as in the 2D case, (^Pe^) , which 

sets the constraint on the choice of Coi analogous to the two-dimensional one (Eq. 
(4.35)). The only significant difference between the 3D and the 2D cases is that in 
3D we cannot demand Co = 0, since the Ei function diverges at zero and the solution 
becomes singular. However, this constraint becomes irrelevant in the “convective” 
regime (Pe' ^ 1), since then c « 1 throughout the bulk of the finger. 

Finally, the finger propagation velocity, derived analogously to the 2D case, is 
given by 


c'(C) 2Pe ^exp(iPe') 

Fmmz = Ei (iPe') - Ei (iPe'Co^)' 


(6.13) 


“^The exponential integral function is defined as Ei(a:) 
^Ei(ai) = 7 + Inai + ai + 0{x'^) 


oo 

— J t~^e~^dt. 

— X 
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with the limiting cases characterized by 

lim i/Pe = —(In^o)”^ (6-14) 

Pe'-i-O 

lim U = K. (6.15) 

Pe'—>cx3 


Again, there is a significant similarity with the 2D case. The (dimensional) propaga¬ 
tion velocity of dissolution fingers at small Peclet numbers is inversely proportional 
to p, 


IVq _ l4>iD 
Peln^o plnCo’ 


(6.16) 


whereas at high Peclet numbers the propagation velocity becomes independent of the 
radius of curvature and yields 


U ^ 7u„. (6.17) 

7. Summary and conclusions. In this paper, we have been analyzing the ge¬ 
ometry of the dissolution fingers, which are commonly found in soluble rocks. We 
have adopted the thin-front approximation, in which the reaction is assumed to take 
place instantaneously with the reactants fully consumed at the dissolution front. In 
this way the problem have become Stefan-like with the undissolved phase downstream 
and the fully dissolved phase upstream. In such a setup we have posed the problem of 
existence and shape of steadily growing fingers. Two potentially limiting assumptions 
have been made. First, we have adopted the Ivantsov ansatz (4.3), which has been 
successfully used in other studies of dendritic growth [28, 24, I], usually leading to the 
parabolic forms of advancing dendrites. Moreover, we have assumed that the trans¬ 
formation (x, y) —>■ (^, r]) from the original Cartesian coordinates to the coordinates 
spanned by the isolines and gradient lines of the reactant concentration field c is con¬ 
formal. Both of these assumptions limit the class of solutions that can be found with 
this technique, leading however to significant simplifications in the analysis, since the 
Ivantsov equation (4.6) is conformally invariant and takes a particularly simple form 
(4.7) in the new coordinate system, where c = c(^). The analysis of these equations 
leads to the conclusion that (^, ij) are parabolic coordinates and thus the advancing 
front is also of a parabolic form. The flow within the finger turns out to be uniform, in 
agreement with the earlier studies on the refraction of groundwater flow on parabolic 
inclusions [31]. The magnitude of the flow is solely the function of the permeability 
ratio between the dissolved and undissolved matrix, n. The concentration field is also 
relatively simple, and can be expressed in terms of the imaginary error function (in 
2d, Eq. 4.31) or the exponential integral (in 3d, Eq. 6.II). This time, the parameter 
controlling the shape of the field lines is Pe^, defined as the ratio of convective to 
diffusive fluxes within the wormhole on the length scale of the radius of curvature of 
the finger tip. 

However, these simplifications come at a price, as we can only impose the con¬ 
centration boundary conditions on the ^ = const lines, which limits the usefulness of 
these solutions in the physical applications. The most relevant physically seems to 
be the case of large Peclet numbers, since in that case the concentration within the 
finger becomes uniform with an exception of thin boundary layer at the reaction front 
itself and the exact placement of Dirichlet boundaries becomes irrelevant. 

Finally, we have obtained the finger propagation velocity, which turns out to be 
inversely proportional to the radius of curvature in the small Pe^ limit and constant for 
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large Pe^ Interestingly, the former prediction seems to be confirmed by the analysis 
of the shapes of the solution pipes — the inspection of Fig. 1.1 and other photographs 
from this site reveals that the longer fingers have invariably smaller radii of curvatures. 

The scale invariance of the problem leads to the dynamical indeterminacy: any 
advancement velocity is permissible provided that the radius of curvature of the finger 
is appropriately tuned, there is no selection principle to fix the tip curvature and the 
finger propagation velocity independently. Finding such a selection principle is in 
general a highly nontrivial task [43] and involves the introduction of additional length 
scale into the system. In the case of viscous fingering and dendritic growth this is 
achieved through the short-scale regularization mechanisms such as surface tension 
or kinetic undercooling. It is not clear what would be the physical origin of such a 
short-scale regularization in the present case. A candidate for the additional length 
scale could be the reaction front width, which in our analysis is zero, but remains 
finite in any physical case. Another possibility is that the additional length scale 
sought here is connected with the presence of other fingers, as it is always the case in 
natural systems (c/. Fig. 1.1). The full solution of the problem would then require 
matching the solution near the parabolic tip with that in the neighborhood of the 
root of the finger, where it joins with its companion. 
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